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QUANTILE  ESTIMATION  IN  DEPENDENT  SEQUENCES 


by 

P.  A.  W.  Lewis 
Naval  Postgraduate  School 
Monterey,  CA  93940 

ABSTRACT  l  '■ 

i 

Standard  nonparametric  estimators  of  quantiles  based  on 
order  statistics  can  be  used  not  only  when  the  data  are  i.i.d., 
but  also  when  the  data  are  from  a  stationary,  ^-mixing  process 
of  continuous  random  variables.  However,  when  the  random  vari¬ 
ables  are  highly  positively  correlated,  sample  sizes  needed  for 
acceptable  precision  in  estimates  of  extreme  quantiles  are  com¬ 
putationally  unmanageable.  A  practical  scheme  is  given,  based 
on  a  maximum  transformation  in  a  two-way  layout  of  the  data, 
which  redqces  the  sample  size  sufficiently  to  allow  an  experi¬ 
menter  to  obtain  a  point  estimate  of  an  extreme  quantile.  Three 
schemes  are  then  given  which  lead  to  confidence  interval  esti¬ 
mates  for  the  quantile.  One  uses  a  spectral  analysis  of  the 
reduced  sample.  The  other  two,  averaged  group  quantiles  and 
nested  group  quantiles,  are  extensions  of  the  method  of  batched 
means  to  quantile  estimation.  None  of  the  schemes  requires  that 
the  process  being  simulated  is  regenerative. 


P.  Heidelberger 

IBM  Research  Center 

York town  Heights,  NY  10598 


1 


9/15/81 


QUANTILE  ESTIMATION  IN  DEPENDENT  SEQUENCES 


1.  INTRODUCTION 

Let  {Xn>  be  a  sequence  of  random  variables  which  is 

assumed  to  be  stationary  with  marginal  distribution 

F  (x)  =  P{X  <  x}  .  For  any  p  ,  0  <  p  <  1  ,  the  p-quantile 

a  n  1 

of  Fx(x)  satisfies  the  equation  p  =  Fx(xp)  if  Fx(x)  1138 
a  positive  density  in  a  neighborhood  of  x  .  This  paper  is 

r 

concerned  with  the  estimation  of  quantiles  xp  from  samples 

Xn  ,  n=l , . . . ,N.  In  particular  the  paper  is  concerned  with  both 

point  estimates  and  confidence  interval  estimates  for  quantiles 

when  the  X  's  are  highly  correlated.  This  is  the  usual  case 
n 

when  the  xn's  are  the  outputs  from  a  simulation  of  a  stochastic 
system,  for  example  the  waiting  times  in  a  queue. 

in  the  case  of  independent  xn's  '  methods  for  quantile 
estimation  are  well  known  (e.g.  Connover,  1980,  pp.  71).  Thus 


let  X,  .  denote  the  order  statistics  from  the  sample  of  size 
(n) 

N  ,  and  let  |_xj  denote  the  greatest  integer  less  than  or  equal 
to  x  .  Then  a  standard  nonparametric  estimator 


Xp  X(|_pN  +  lj) 

is  a  consistent,  asymptotically  normal,  estimator  of  x  with 
bias  which  is  0(N-1)  and  variance  [p (1-p) ]/[ f2 (xp) N]  +  0(N  2)  , 
where  fv(x  )  is  the  derivative  of  Fv(x)  at  x  .  This  vari- 

X  p  A  p 

ance  can  be  estimated  by  estimating  fv (x  )  ,  but  nonparametric 

A  P 

confidence  intervals  can  also  be  obtained  for  x^  using  the 
order  statistics  of  the  sample  (Connover,  1980,  pp.  111-116). 
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The  difficulty  with  these  estimation  methods  is  that 

they  require  a  large  amount  of  computer  storage  and  computing 

time  to  sort  the  sample.  Also  for  high  or  low  p  ,  the  order- 

statistic  estimator  might  be  biased  and  very  non-normal.  A 

solution  is  to  use  the  maximum  transformation  (Goodman,  Lewis 

and  Robbins,  1971).  This  not  only  moves  the  problem  back  to  one 

of  estimating  a  median,  but  also  allows  the  use  of  stochastic 

approximation  (Robbins -Monro)  methods.  A  very  satisfactory 

method  of  this  type  using  a  minimal  amount  of  memory  and  giving 

both  point  and  confidence  interval  estimates  for  the  unknown 

quantile  x  has  been  given  by  Robinson  (1975). 

P 

For  the  case  of  dependent  *n's  •  the  usual  situation 
encountered  in  system  simulation  studies,  quantile  estimation 
is  much  more  difficult  than  in  the  independent  case.  The  point 

A 

estimate  x  is  still  valid;  however  its  variance  is  inflated 
P 

by  a  factor  p(0;x  )  .  Here  p(0;x  )  is  the  initial  point  on 

P  P 

the  spectrum  of  the  binary  process  {I  (x  ) }  ,  defined  for  each 

**  Jr 

n  to  be  1  if  X  <  x  and  zero  otherwise.  More  directly  it  is 

n  -  p 

P(0 ; x  )  =  lim  n  Var{(I  (x  )  +..+  I  (x  ))/n)  . 

P  n-*-°°  "  P 


The  nonparametric  confidence  interval  estimation  procedures  are 
no  longer  directly  applicable.  However  it  is  possible  to  esti¬ 
mate  p(0;Xp)  using  methods  of  Heidelberger  and  Welch  (1980, 

1981)  and  to  estimate  the  density  f (x  )  by  standard  methods; 

P 

used  together  these  give  an  estimate  of  the  variance  of  xp 
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and  large-sample  confidence  intervals.  The  main  problem  with 
these  order-statistic  estimates  is ,  however,  that  the  sample 
sizes  required  to  obtain  reasonable  precision  with  the  estimates 
are  prohibitive  when  the  Xn ' s  are  highly  positively  correlated . 

The  basic  scheme  considered  in  this  paper  to  handle  the 
sample  size  problem  for  dependent,  <(>-mixing  sequences  is  to  set 
out  the  data  in  a  (conceptual)  v  *  m  array 


Xk,i  Xi  +  (k-l)m  ' 


for  k  =  l,...,v  and  i  =  1,...  m  =  N/v  ,  where  v  is  often 
but  not  necessarily  chosen  so  that  pv  =  0.5  .  Then  if  m(and  N) 
is  large  enough  so  that  Xn’s  which  are  m  apart  are  approxi¬ 
mately  independent,  the  table  can  be  collapsed  by  taking  maxima 
down  the  columns.  The  maxima 


Y.  =  max  X.  .  , 

1  l<k  <v  K,;L 


for  i  =  l,...,m  are  now  a  reduced  sample  whose  q=p  quantile,  y 
corresponds  to  the  desired  quantile,  Xp  ,  of  the  {Xfi}  process. 
The  main  point  of  this  procedure  is  that  it  gives  a  very  sub¬ 
stantial  sample  size  reduction.  In  addition,  there  is  generally 
slightly  less  overall  correlation  in  the  Y^  sequence  and  this 
counteracts  the  increase  in  variance  which  occurs  with  the  use 
of  the  maximum  transformation  and  order  statistic  estimates  of 
the  quantile  x  . 


Since  the  Y^'s  (the  maximum  transformed  sample)  are 
still  correlated  several  problems  remain,  notably  estimating 
the  variance  of  the  point  estimate  'the  q=pV  sample  quantile  of 
the  Y.'s)  and  obtaining  confidence  interval  estimates  of  x  . 

X  H 

Three  confidence  interval  schemes  are  considered;  averaged  group 
quantiles,  nested  group  quantiles  and  a  scheme  based  on  esti¬ 
mating  the  probability  density  function  of  the  Yi's  at  xp  , 
and  the  initial  spectral  point,  p(0;yg)  of  the  binary  process 
obtained  from  comparison  of  the  Y^'s  to  ,  i.e.  I^(yg)  . 

Extensive  empirical  sampling  studies  using  M/G/l  queues 
and  stationary  sequences  of  autocorrelated  exponential  random 
variables  show  that  the  above  schemes  provide  a  reliable  method 
for  estimating  a  quantile  in  a  dependent  sequence  { > 


2.  QUANTILE  ESTIMATION  AND  THE  MAXIMUM  TRANSFORMATION 
1 .  Quantile  estimation  for  i.i.d.  sequences. 

Let  X1 , .  .  .  ,Xn,..  .  ,XN  be  a  sample  of  i.i.d.  random 

variables  from  a  continuous  distribution  Fx(x)  with 
probability  density  function  F  (x)  •  For  0  <  p  <  1  let 


(2.1)  xp  =  inf{x  :  Fx(*)  =  P^  =  Fx  ^  ' 


where  F^fp)  is  the  inverse  of  Fx(x)  with  derivative 

1/f  (x  )  .  The  quantity  x  is  the  pth  quantile  of  Fx<*) 

A  p  H 


<  +a»  jr*  >  *rv**v.«*.»  r* 


V. 


with  L  <  U  chosen  in  such  a  way  as  to  make  the  probability  as 
close  as  possible  to  the  desired  coverage  (1-a)  . 
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Another  method  for  estimating  quantiles  uses  the  method 
of  stochastic  approximation  (Robbins  and  Monro  (1981)).  This 


provides  a  method  of  quantile  estimation  which  avoids  sorting 
and  requires  only  a  minimal  amount  of  storage.  The  asymptotic 
variance  of  the  stochastic  approximation  estimate  is  given  by 
(2.4).  A  very  satisfactory  development  of  this  method  is  given 
in  Robinson  (1975) .  Since  it  is  not  extendable  to  dependent 
data,  it  is  not  described.  The  key  idea  in  the  method,  however, 
is  the  use  of  the  maximum  transformation,  and  since  this  is  of 
use  in  the  present  context,  it  is  described  next. 

2.2.  The  maximum  transformation 

The  maximum  (minimum)  transformation  is  a  computationally 
simple  transformation  and  compaction  of  i.i.d.  data  which  trans¬ 
forms  the  problem  of  estimating  an  extreme  quantile  of  a  random 
variable  into  that  of  estimating  a  more  central  quantile,  e.g. 
the  median.  Thus  assume  that  p  >  1/2  ,  and  let 

Y  =  max  (X-^ ,  .  .  .  ,Xv)  .  (When  p  <  1/2  a  minimum  transformation  is  used). 
Then 

(2.7)  P{Y  <  x  }  =  F  (x  )  =  (F  (x  ) }V  =  pV  =  q  . 

F  1  P  A  P 

Thus  the  pth  quantile  of  Fx  is  the  q  =  pvth  quantile  of  the 

Y  variable.  In  particular  let  v  =  (_£n  (1/2)  /£n  (p)  J  ;  then  x 

P 

is  approximately  the  median  of  the  Y  variable.  Details  are 
given  in  Goodman,  Lewis,  and  Robbins  (1971). 

There  is  a  price  paid  with  this  transformation  in  infla¬ 
tion  of  the  leading  term  of  the  variance  (2.4).  Thus  let  v 
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divide  N  ,  so  that  N/v  —  m  .  Then  the  X  sample  is  reduced 
to  a  sample  by  application  of  the  maximum  transfor¬ 
mation  to  each  successive  non-overlapping  gxoup  of  X, 's  of 


size  v  .  Then  the  density  of  the  Y/s  is  f^(xp)  =  v  f^(x  )p 
and  the  variance  of  the  order  statistic  estimator, 


v—  1 


q  ' 


of  x_  in  the  sample  Y,  , . .  .  ,Y  is 
P  i  m 


(2.8) 


Pll-P) 


var(y  )  - 

NfX(V 


(1-PV) 


v  (1-p)  p 


v—  1 


var(x  ) 


(1-PV) 


v(l-p) p 


V—  1 


The  right-hand  multiplicative  factor  in  this  last  expression  is 
approximately  1.4  if  v  is  chosen  to  make  q  =  pv  approximately 
0.5  . 

When  using  stochastic  approximation,  shifting  to  the  neigh¬ 
borhood  of  the  median  is  essential  for  the  method  to  be  a  well- 
behaved  estimation  procedure.  For  estimating  several  quantiles 
simultaneously,  the  method  is  applicable  in  a  nested  scheme 
which  is  very  efficient  with  respect  to  storage  and  speed. 

For  p  <  1/2  a  minimum  transformation  is  used.  Next- 
to-maximum  and  maximum-minimum  schemes  give  more  flexible  and 
robust  schemes  but  are  not  discussed  in  the  present  paper. 

2.3.  Dependent  Data;  General  Considerations  and  Order  Statistics 

Quantile  estimation  in  dependent  data  is  an  order  of 
magnitude  more  difficult  than  in  independent  data.  Fishman 
(1978,  p.  270  )  notes  that  no  satisfactory  solution  exists. 
Iglehart  (1976),  Seila  (1976)  and  Moore  (1979)  have  given  special 
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=  P(Xn  <  X  ,  Xntk  <  x}  -  P(Xn  <  x)  P{xn+k  <  X)  . 

The  notation  p(0;xp)  comes  from  the  fact  that  this 

quantity  is  the  initial  point  (f=0)  of  the  spectrum  of  the 

process  {I  (x  ) }  ; 

n  p 

OO 

(2.13)  p(f;x)  =  l  cos(2-rrfk)  y,  (x)  ,  -  i  <  f  <  1/2 

k=-°°  K  z  ~ 

Now  the  Central  Limit  Theorem  for  x  =  X , ,  . i , 

p  ( |  np+lj  ) 

(Sen,  1972)  states  that  if  fx(x)  continuous/ 

bounded  and  non-zero  in  some  small  neighborhood  of  x  ,  and  if 

I 

f^(x)  is  bounded  in  this  neighborhood,  then 

.,1/2  .  ~ 

N  '  (x  -  x  ) 

(2.14)  - £  -IP - -  ( 0 , 1) 

(p(0;xp))1/^/fx(xp) 

if  the  process  { Xn >  is  <})-mixing  and 

(2.15)  ol  =  p (0 ;x  ) /{fx(x) }2 

Ap  P  A  P 

is  finite.  For  details  on  <f-mixing  see  Sen  (1972)  and 
Billingsley  (1968);  some  discussion  is  given  in  Section  2.4. 

Regenerative  processes  such  as  the  M/G/l  queuing 
system  waiting  time  (wn^  are  $~nuxing. 

The  problem  with  using  xp  as  an  estimator  in  positively 
correlated  sequences  is  that  the  sample  sizes  required  for 
adequate  precision  are  prohibitively  large.  Both  sorting  times 


V, 


and  memory  times  are  then  unrealistic.  A  measure  of  the  infla¬ 
tion  of  sample  size  over  the  independence  case  is  the  ratio  of 
p ( 0 ; Xp )  to  its  value  p(l-p)  for  the  i.i.d.  process  with 
identical  marginal  distributions.  This  has  been  investigated 
by  Blomqvist  (1967)  for  the  M/G/l  queue.  For  extreme  quantiles 
of  the  M/M/1  queue  with  traffic  intensity  p  =  0.9  this  ratio 
is  400  .  Greater  ratios  are  possible,  depending  on  the  traffic 
intensity  and  the  skewness  of  the  service  time  distribution. 
Specifically,  the  M/M/1  queue  with  p  =  0.9  requires  a  sample 
size  of  roughly  500,000  customers  to  estimate  the  0.99  quantile 
of  the  waiting  time  distribution  to  within  plus  or  minus  10% 
accuracy.  This  is  derived  from  Table  8  of  Blomqvist  (1967) ;  more 
precisely  this  is  the  sample  size  required  for  a  90%  confidence 
interval  for  x  to  have  a  relative  half-width  of  0.10  .  For 
the  0.999  quantile  the  required  sample  size  is  approximately 
2,300,000.  Clearly  storing  and  sorting  the  entire  sequence  is 
impractical  in  such  cases.  Actually  to  produce  an  order-statistic 
point  estimate  of  x  requires  storing  only  the  largest  (l-p)N 

r 

values  of  the  sequence.  However,  this  ordering  must  be  dynam¬ 
ically  maintained  as  the  sequence  is  generated,  a  computationally 
expensive  operation.  Additional  storage  is  required  to  estimate 

the  variance  term  p(0;x  ):  the  positions  in  the  sequence  at 

P 

which  the  (l-p)N  largest  values  occur  must  also  be  saved  in 
order  to  construct  the  binary  sequence  {I^fx^)}  •  Furthermore, 
the  Discrete  Fourier  Transform  of  the  extremely  long  sequence 
{In(x  )  ,  n  =  1,...,N)  must  be  taken  to  accomplish  the  variance 

estimation . 
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Some  of  the  storage  and  sorting  problems  could  be  re¬ 
lieved  by  using  stochastic  approximation  .  However  properties 
of  the  stochastic  approximation,  particularly  its  asymptotic 
variance,  are  unknown  for  dependent  data.  More  importantly 
direct  application  of  the  maximum  transformation  to  bring  the 
quantile  estimation  down  to  a  problem  of  estimating  the  median 
requires  independence  of  the  xn ' s  • 

We  now  present  a  method  for  using  the  maximum  transfor¬ 
mation  t.  achieve  sample  size  reduction  and  a  practical  scheme 
for  point  and  confidence  interval  estimates  with  dependent  data. 

2.4.  The  maximum  transformation  in  the  dependent  case 

The  basic  idea  behind  the  use  of  the  maximum  transfor¬ 
mation  is  to  combine  elements  of  X  ,  n  =  1 ,2  , .  .  .  ,  in  a  <j>- 

mixing  process  which  are  sufficiently  far  apart  so  as  to  be 

approximately  independent.  To  define  <|>-mixing  let  and 

Mn+m  be  respectively  the  o-fields  generated  by  {Xi;  i  <  n} 

and  (X.  ;  i  >  n  +  m}  .  If  E.eMn  and  E_eM°°.  ,  then 

x  —  i  — °°  2.  n+m 

{ Xn }  is  Q-mixing  if  for  all  n(-«  <  n  <  ®)  and  m(>^  1) 

(2.16)  |  P  (E2 1 E1)  -  P(E2)|<  <Mm)  ,  <f>  (m)  >  0  , 

where  1  >  <J>(1)  >  <f>(2)  >_.  ..  and  lim^^  $  (m)  =0  . 

Thus  we  set  out  the  data  in  a  v*m  array,  where 

m  =  N/v,  as 
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(2.17) 


Xk,i  Xi+(k-l)m 


k  1 » • • . / v ;  i-1 ,  .  .  .  ,m 


We  assume  for  the  moment  that  elements  in  the  m  columns 
are  independent;  since  we  are  assuming  that  v  is  fixed  and 
the  process  (Xn)  is  4>-mixing  this  will  certainly  be  true  as 
N  (and  m)  get  very  large.  Now  applying  the  maximum  transform 
down  the  columns  we  get  a  reduced  series  Y^  ,  where 


(2.18)  Y.  =  max  X,  .  i  =  1,2,. . . ,m  . 

1  l<k<v  K'1 


Values  of  v  required  to  reduce  the  problem  to  that  of  esti¬ 
mating  the  median  of  the  Y^s  are  given  in  Table  2.1  . 


Size  of 

Table  2. 

v  for 

1 

PV  =  1/2 

p 

0.900 

0.950 

0.975 

0.980 

0.990 

0.995 

0.998 

0.999 

V 

6.6 

13.5 

27.4 

34.3 

69.0 

138.3 

346.2 

692.8 

We  emphasize  here  that  transformations  to  points  other 

than  the  median  are  possible,  and  are  essential  when  considering 

simultaneous  estimation  of  several  quantiles  or  estimation  of 

the  median  of  the  X  's  .  However  in  what  follows  transforma- 

n 

tion  to  the  median  is  usually  assumed.  There  will  be  an 
attendant  bias  because  pv  will  never  quite  be  1/2  . 

The  purposes  of  the  maximum  transformation  are 
(i)  to  reduce  the  attendant  sample  size.  It  should  be  kept  in  mind 


13 


that  for  higher  quantiles  larger  sample  sizes  are  needed  for 
given  precision.  The  result  is  that  for  fixed  precision  the 
series  is  roughly  the  same  length  for  all  p  >  0  .  In 
particular  for  p  =  0.999  the  reduction  of  the  sample  size  is 
by  a  factor  of  693  which  is  sufficient  to  reduce  a  series  of 
completely  unmanageable  length  to  one  which  can  clearly  be 
accommodated  in  a  digital  computer  (see  Section  5) . 

(ii)  To  reduce  the  problem  to  estimation  of  a  median  rather 
than  an  extreme  quantile.  Since  stochastic  approximation  is 
not  used  with  the  Y^  series,  this  is  not  essential  for  point 
estimates.  It  is,  however,  helpful  in  obtaining  confidence 
interval  estimates. 

(iii)  To  possibly  reduce  correlation.  It  will  be  shown  that 
the  correlations  in  the  {Y^}  series  are  slightly  less  than 
the  correlations  in  the  {Xn>  series.  Though  a  small  effect, 

it  is  usually  sufficient  to  offset  inflation  in  variances  due  to 

use  of  the  maximum  transformation.  To  see  this  assume  that  the 

structure  of  I  (x  )  is  Markovian  with  p(0;x  )  ~  (l+p)/(l-p)  . 

n  p  p 

If  p  =  0.99  this  equals  199.  If  p  for  the  Y^  process 
chopped  at  its  median,  i.e.  {ln(xQ  ^ ) }  ,  is  reduced  to  only  0.98, 
then  the  new  p(0;x  )  ~  99  . 

r  p  ~ 

Clearly  if  m  is  too  small,  the  assumption  of  inde¬ 
pendence  down  the  column  will  be  violated.  Some  analysis  of 
this  effect  is  possible.  Expanding  the  definition  (2.18)  to 
account  explicitly  for  the  parameter  v  ,  the  number  of  Xn ' s  , 
v  apart,  over  which  the  maximum  is  taken,  so  that 

Y.(v)  =  max  X,  .  ,  i=l,2,...,m  ,  we  have  from  Billingsley 

1  l<k<v  K>1 

(1968,  pp.  174,  20.49)  that 
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(2.19) 


|P{Yi(v)  <  xp}  -  pV}|<  v0(m)  . 


Thus  we  want  m  as  large  as  possible  and,  for  given  v  ,  the 
difference  in  the  probabilities  goes  to  zero  as  m  -+  »  ,  by  the 
definition  of  mixing.  A  slightly  tighter  bound  is  given  by  the 

v-2  , 

Lemma.  [p{Y.;  (v)  <  x^l  -  pV)  <  (m)  p  j  p 

p  k=0 

for  v  >_  2  . 

A  straight  forward  proof  of  this  is  obtained  using  induction. 

The  point  estimate  used  for  x^  is 

A 

yq  ~  Y  ( [_mq+lj  +1) 

where  q  =  pV  .  Since  the  Y^'s  are  dependent  the  properties 
of  this  estimator  are  given  as  at  (2.14)  and  (2.15)  ,  with 
proviso  for  some  peculiarities  in  the  structure  of  the 
process  which  we  discuss  now.  Finite  sample  properties  of  the 
estimator  y  are  studied  by  simulation  in  Section  5. 

2.5.  The  maximum- transformed  process  (Y^ }  . 

The  maximum- transformed  process  •  obtained  from 

the  process  by  taking  the  maximum  of  xn's  which  are  m 

apart,  is  not  a  stationary  process  in  the  same  sense  that  the 

{X  }  process  is  stationary.  To  see  this  consider 
n 
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(2.20) 


Y.  =  max  X,  .  =  max{X. ,X  .. 
1  l<k<v  kfl  1  m+1 


'X(v-l)m+l}  ' 


It  will  be  correlated  with  successive  Y^'s  *  and  this  correla¬ 
tion  can  be  expected  to  decrease  since  the  {Xn>  process  is 
mixing.  However,  the  correlation  will  eventually  increase  because 

(2. 2D  Y  =  max{X  ,X_  , . . . ,X  } 

m  m  2m  vm 

is  coupled  to  Y^  by  the  v-1  dependent  pairs 

{Xm'Xm+l^ '*•*  {X(v-l)m'X (v-1) m+l}  *  Thus  the  {Y.}  process 
is  circular.  However,  it  is  not  stationary  in  a  circular  sense 
because  the  correlation  between  Y^  and 

Y2  =  max{X2 ,Xm+2 , . . . ,X m+2 1  is  through  the  v  dependent 

pairs  {X1,X2},  txm+1'xm+2} ' • • • '{X(v-i)m+l'X(v-l)m+2}  '  and  is 

therefore  different  from  the  correlation  between  Y.  and  Y 

l  m 

In  a  circularly  stationary  process  these  correlations  would  be 

the  same.  Towards  the  middle  of  the  process,  say  i  =  m/2  , 

the  lag  one  correlations  will  be  the  same  and  there  will  be  no 

"edge  effect"  because,  if  m  is  large  enough,  the  correlation 

between  Y  and  Y.  ,  and  Y  and  Y  will  be  zero, 

m/ 2  1  m/ 2  m 

Without  belabouring  this  "edge  effect"  and  circularity 
in  the  { Y^ }  process  there  are  two  main  points  to  be  made. 

(i)  For  finite  m  it  will  introduce  bias  because,  for 
instance,  methods  using  spectral  techniques  for  estimating 
p(0;x)  in  Section  3  are  based  on  assumptions  of  stationary. 


Note,  however,  that  {Y^}  is  marginally  stationary,  so  density 
estimation  techniques  for  f^(x)  are  not  affected. 

(ii)  For  m  large  the  effect  will  be  negligible  and 
the  asymptotics  go  through  in  the  usual  way  by  ignoring  a  fixed 
(or  slowly  increasing)  set  of  Y^'s  at  the  beginning  and  end  of 
the  process.  The  usual  assumption  is  that  the  extent  of  this 
set  is  smaller  than  the  extent  of  the  dependence  in  the  (Yi) 
process. 
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_ _ 

y 

3.  CONFIDENCE  INTERVALS  FOR  QUANTILES  IN  DEPENDENT  DATA 

In  the  previous  section  we  gave  a  point  estimator  y 

q  I 

for  the  quantile  of  the  marginal  distribution  in  a  stationary  ' 

time  series.  In  this  section  we  consider  three  methods  for  t  1 

generating  confidence  interval  estimates  for  these  quantiles. 

The  first  method  is  an  extension  of  the  spectral  method  of 
Heidelberger  and  Welch  (1980;  1981).  The  second  and  third  methods 
are  extensions  of  the  method  of  batch  means  (see  e.g.  Law  and 

I 

Carson  (1979)  or  Fishman  (1978)).  These  last  two  methods  also 

! 

give  two  new  point  estimates.  j 

These  methods  for  generating  confidence  interval  esti¬ 
mates  do  not  rely  upon  use  of  the  max-transform  and  the  case 
v  =  1  corresponds  to  working  with  the  original  time  series. 

Moreover,  if  the  max-transform  is  used  we  do  not  require  that  j 

V  *  * 

p  =  0.50,  i.e.  v~  In. 5/ In  p  .  There  are  considerable  compu¬ 
tational  savings,  however,  from  using  some  maximum  transformation. 

j 

This  and  statistical  considerations  will  generally  dictate  a  j 

maximum  transformation  using  a  v  slightly  less  than  the  median 
transform  v  .  Another  determing  factor  is  that  one  generally 
wants  to  estimate  more  than  one  quantile.  This  factor  will  be 
considered  elsewhere. 

3.1.  A  Spectral  Method 

Let  { Y ^ , . . . , Ym }  be  the  sequence  of  max-transformed 
variables  as  defined  by  equation  (2.18).  The  point  estimate  of 
Xp  is  the  qth  order  statistic  of  {Y^,...,Ym}  '  Yq '  w^ere 
q  =  pv  .  Let  f^(x)  be  the  stationary  density  function  of  Y^ 
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V, 


and  let  p(0;x)  be  the  stationary  spectral  density  of  the  in¬ 
dicator  function  process  I^(x)  associated  with  {y^}  at  zero 
frequency  as  defined  by  equation  (2.13).  Here  edge  effects 
because  of  the  quasi-circularity  of  the  Y^  process  are  ignored. 
For  large  values  of  m  (see  Sen  (1972)  and  Babu  and  Singh  (1978)) 
/in  (y  -  x  )/(p(0;x  )/f  (x  ))  '  has  a  normal  distribution  with 

Mr  Ir  y  r 

mean  zero  and  variance  one.  Confidence  intervals  for  x  based 

P 

on  this  Central  Limit  Theorem  are  generated  by  estimating  both 

p(0;x  )  and  f  (x  )  at  the  estimated  quantile  .  This  is  the 
p  y  p 

first  method  alluded  to  above. 

The  quantity  p(0;x  )  is  estimated  using  the  spectral 

P 

method  of  Heidelberger  and  Welch  (1980;  1981)  applied  to  the 
sequence  (i^(yq)  ,  i=l,...,m}  .  This  method  uses  least  squares 

to  fit  a  low  order  polynomial  of  degree  d  to  the  logarithm  of 
the  first  K  values  of  the  averaged  periodogram  of  d^(y  )}  . 

As  suggested  in  Heidelberger  and  Welch  (1980;  1981)  we  used 
K  =  25  and  d  =  2  ,  although  for  extreme  quantiles  in  highly 
congested  queues  with  short  run  lengths,  d  =  3  was  required 
to  produce  valid  confidence  intervals.  This  point  is  discussed 
more  fully  in  Section  4. 

We  found  little  or  no  loss  in  confidence  interval  cov- 

A 

erage  by  using  the  estimated  quantile  y  with  the  sequence 

as 

(I^(yq)  }  rather  than  using  the  known  quantile  x^  with  the 
sequence  (I^(Xp))  .  However,  a  proof  of  the  convergence  of 
the  distributional  properties  of  the  periodogram  of  (I^fy^)} 
appears  difficult. 
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The  density  f^(Xp)  t^ie  niax-transformed  variables 

Yl'’’*,Ym  may  be  estimated  using  a  standard  kernel  density 
estimator  (see  Parzen  (1962)  and  Rosenblatt  (1956)).  Specifically, 
let  W(x)  be  a  weighting  function  such  that 


0  <  W  (x)  <  W  <  “  , 


(3.1) 


/  W(x)dx  =  1  , 


lim  | xW (x) |=0. 
x  I  ■>«> 


Let  b(m)  be  a  sequence  of  bandwidth  constants  satisfying 
lim  b(m)  =0  and  lim  mb(m)  =  00  .  For  any  x  the  kernel 

m-*oo  m->°° 

density  estimate  of  fy(x)  is  defined  by 


(3.2) 


^  -j  m  1 

f„(x)  =  ±  l  W((x-Y.)/b(m)  )  . 


m  b (in) 


1/2 

Under  the  conditions  (3.1)  and  if  l  (n )  '  <  00  ,  where  the 

n=0 


<j>(n)'s  are  defined  at  (2.16)  ,  then  fy(x)  converges  in  prob¬ 
ability  to  fy(x)  . 

In  our  case  we  require  a  density  estimate  at  the  unknown 

A  A 

point  x^  .  Simple  Taylor  series  expansions  show  that  fv(y„) 
-  p  y  4 

converges  in  probability  to  fy(x^)  if  the  first  k  derivatives, 


W 


(k) 


(x)  ,  of  W(x)  satisfy 


/  |  W  ^  ^  (x)  |  dx  <  00  and  lim  |xw'*;  (x)  |  =  0  and  if 

—00  I  x  I  •♦■oo 


W{k)  (x)  |  <  W  <  00  ,  if 
,  (k) 
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b(m)  =  Zmm  where  Zm  converges  in  probability  to  a  positive 
constant  and  0  <  c  <  k/2(k+l)  .  If  the  weighting  function 

W(x)  does  not  satisfy  the  different  ability  conditions  then 
more  delicate  arguments  are  required  to  show  convergence  of 

A  A 

fy(Xp)  (see,  for  example,  Robinson  (1975)). 

We  have  experimented  with  two  weighting  functions,  a 
triangular  window 


(3.3) 


and  a  cosine  window 


x  |  £  1  , 
x|  >  1  , 


(3.4) 


W(x) 


j 1/2  cos(x) 

to 


I  x  |  <  tt/2  , 

|  X  |  >  tt/2  ; 


Bandwidth  sequences  b (m)  had  the  form  b(m)  =  Z  m 

m 

for  various  powers  of  0  <  c  <  1/2  and  random  variables  Z 

m 

which  ranged  from  constants  to  measures  of  the  spread  of  the 

distribution  of  the  Y^'s  .  We  found  the  density  estimates  to 

be  relatively  insensitive  to  both  the  shape  of  the  weighting 

function  and  the  parameters  of  the  bandwidth  sequence.  For 

small  samples  they  tended  to  slightly  overestimate  fy(x^)  but 

converged  to  fy(x^)  for  large  samples.  In  Section  4  we  report 

the  results  of  experiments  using  the  triangular  window  with 

c  =  1/3  and  Zm  equal  to  the  inter-quantile  range  of  the 

Y.'s  ,  i.e.  b  (m)  =  <Y,,  ..  -  Y,,  ,.l  m  . 

J  X  ( L . 75m+l J )  ( L . 2  5m+l J ) ( 
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3.2.  Averaged  and  Nested  Group  Quantiles;  Definition 


We  now  describe  two  methods  of  generating  confidence 
intervals  for  the  quantile  x^  which  do  not  require  direct  es¬ 
timation  of  p(0;Xp)  for  the  {Y^}  sequence.  These  methods, 
which  we  call  averaged  group  quantiles  (agq)  and  nested  group 
quantiles  (ngq) ,  are  extensions  of  the  method  of  batch  means  to 
quantile  estimation.  Seila  (1976)  considered  a  version  of  agq 
without  the  max  transform  (v=l)  for  regenerative  processes. 

Let  N  be  the  total  sample  size  and  divide  {X^,...,XNJ 
into  G  non-overlapping  groups  with  mv  observations  in  each 
group  (N  =  Gmv)  .  Define 

(3'5)  X£,j,k  =  X(£-l)mk+j+(k-l)m  for  £=1,...,G, 

j=l, • • • ,m, 
k=l , . . . , v. 


The  subscript  £  refers  to  the  group  number,  and  the  data  in 
each  group  can  be  thought  of  as  forming  a  matrix  of  dimension 
v  by  m  .  The  first  row  of  this  matrix  is  formed  from  the 
first  m  observations  in  the  group,  the  second  row  is  made  up 
of  the  second  m  observations  in  the  group,  etc.  For  any  fixed 
values  of  £  and  j  the  points  X£  j  ^  and  j  k+y  d£k<v) 

are  at  lag  m  (i.e.,  separated  by  m-1  observations)  and  for 
large  values  of  m,  X0  .  . ,  X„  .  ,,..., X  .  ,  are  approximately 

a  ,  J  f  1  x.  ,  J  t  & 

independent.  Define  Y£j  =  max  X£^k  .  thus  y  is  the 
jth  max-transformed  variable  in  the  £th  group.  For  large  m 


i  xp^ ::  q  =  p  with  lp{Yi;,j 

Section  2. 


x  }  -  q|  <  (fi(m)v  as  in 
P 


F°r  eaCh  *=1 . G  let  Y£(l)  i  Y£(2)---I  Y?(m) 

the  order  statistics  of  (Y^  Y  .  ..  Ypmf  dnd  define 

yq£  =  Y£  (  [_mq+lj  )  '  i,e-  yq£  is  an  order  statistic  estimate  of 

xp  derived  from  the  £th  group.  We  call  a  group  quan¬ 

tile  estimate.  As  in  the  case  of  batch  means,  for  large  values 
of  m  ,  { yqi , . . . ,yqG)  wil1  be  distributed  as  i.i.d.  normals 

with  mean  xp  and  variance  p (0 ; xp) /mf ^ (x  )  .  This  suggests 

several  point  and  interval  estimates  for  x 

P 

The  averaged  group  quantile  estimate  is  defined  by 


(3.6) 


1  ^ 
agq(xn)  =  -  l  y 


A  confidence  interval  for  xp  is  formed  by  assuming  that 


(3.7) 


/G  (agq(xp)  -  xp)/S(agq) 


has  a  students-t  distribution  with  G-l  degrees  of  freedom 


where 


1  ^  ^  2 
(agq)  =  azj  l  (ynf  -  agq(x_))  • 


This  confidence  interval  is  valid  if  {y  n  ...,y  }  are  i.i.d. 

ql  t  cjG 

normals  with  mean  x  and  finite  variance. 

P 

The  nested  group  quantile  point  estimate  of  x  is 

defined  to  be  the  median  of  y  ,  ...,y  „  ,  i.e.  if 

ql,  J qG 

<  ...  <■  y.  ate  the  order  statistics  of  v . v  _  , 


then  the  ngm  is  defined  by  (assuming  G  is  odd) 


(3.8) 


ngq(xp)  =  yq(L.5G+1J)  • 


An  approximate  100x(l-a)%  confidence  interval  for  xp  from 
(2.6)  is 


(3.9) 


(Yq(L)  ,  yq(U)}  ' 


where  for  any  L  <  U,  a  is  defined  as  at  (2.6)  ,  by 


(3.10) 


U-l  k 

1  -  a  =  l  (p  (.5)  (.5) 

k=L 


G-k 


This  confidence  interval  is  valid  if  {y  j  • • • /yqG)  are  mutually 

independent  with  P{y  ^  £  xp}  =0.5  .  For  large  samples  this 

is  guaranteed  by  the  asymptotic  normality  of  order  statistics. 

As  m-*°°,  agq(x  )  and  ngqCx  )  are  asymptotically  un- 
P  P 

biased  and  their  confidence  intervals  are  asymptotically  valid, 
since  the  sequence  {X  }  is  assumed  to  be  mixing.  Notice  that 

/v  A 

if  G  =  1  ,  then  agq(xp)  =  ngq(xp)  =  yg  where  yg  is  defined 
by  equation  (2.20),  though  clearly  G  =  1  is  not  reasonable  if 
a  confidence  interval  estimate  is  required.  Furthermore,  as  with 
y^  ,  these  point  and  interval  estimates  are  valid  if  the  max- 
transform  is  not  used  (v=l)  in  which  case  one  is  working  with 
group  quantiles  of  the  original  sequence  {XjP  •  If  the  max- 

v 

transform  is  used  there  is  no  requirement  that  p  = 


0.5  . 


3.3.  Averaged  and  Nested  Group  Quantiles:  Discussion 

There  are  a  number  of  possible  sources  of  error  in  the 

above  two  confidence  interval  procedures.  First,  if  m  is 

small  then,  because  of  dependence,  PiY^  £  xp}  will  differ 

from  q  =  pv  by  a  significant  amount  and  E(y  „)  ^  x  .  Notice 

q  ~  P 

that  with  G  >  1  ,  m  for  ngq(x  )  and  agq(x  )  is  smaller 

Jr  ir 

than  the  m  for  y  and  therefore  this  source  of  error  is 

si 

more  likely  with  ngq(x  )  and  agq(x  )  than  with  y  .  Second, 

p  p  q 

even  if  the  y  ^  sequence  consists  of  i.i.d.  observations,  if 
m  is  small  E(y^)  ^  Xp  due  to  the  bias  in  order  statistics  for 
small  samples.  This  bias  will  be  more  severe  at  the  tails  of  the 
distribution.  This  problem  is  again  potentially  more  serious  with 
ngq (Xp)  and  agq (xp)  than  with  yg  .  Third,  the  random  vari¬ 
ables  y  !  . . . 'YqG  maY  t>e  correlated.  All  three  of  these  fac¬ 
tors  suggest  that  the  number  of  groups,  G  ,  should  be  as  small 
as  possible.  The  choice  G  =  5  is  the  smallest  value  of  G 

^  A 

for  which  P{y  <  y  ,_>  }  >  0.90  (for  L=1,U  =  G  =  5 

q  \ ;  —  p  q  j  — 

the  confidence  level  is  0.9375).  Given  that  G  should  be  small, 

A  A 

any  independence  tests  for  y  ^  •'*,irqG  have  very  low 

power.  The  sampling  experiments  described  in  Section  4  use 
G  =  5  for  both  agq(x  )  and  ngq(x  )  and  the  correlations 

A 

measured  between  the  Vq j ' s  have  been  very  low  even  for  strongly 
correlated  sequences  of  xn ' s  (for  these  experiments  the  sample 
sizes  were  such  that  m  _>  200).  Any  residual  correlations 

between  the  group  medians  could  be  further  reduced  by  deleting 
observations  between  groups,  i.e.  by  discarding  observations 
before  starting  the  next  group. 
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3.3.  Theoretical  Comparisons 

It  is  interesting  to  compare  the  expected  confidence 

interval  widths  using  these  three  methods  under  the  assumption 

that  each  method  is  valid.  Let  m  =  N/v  and  m^,  =  N/vG  =  m/G 

2 

and  assume  (without  loss  of  generality)  that  p(0;x  )/m_f  (x  )  =  1 

p  p 

A  A 

Let  G  =  5  and  assume  that  y  ,  .  ,.,y  _  are  i.i.d.  normals 

1 ql >  qG 

2 

with  mean  x^  and  variance  p(0;x  )/m_f  (x  )  =  1  .  Then  the 
P  P  ^  P 

expected  ngq  confidence  interval  half-width  is  1.163 
(Pearson  and  Hartley  (1966),  Table  28,  p.  190).  Consider  now 
the  agq  method.  The  t-multiplier  for  a  93.75%  confidence  inter¬ 
val  with  4  degrees  of  freedom  is  2.562  so,  assuming  an  unbiased 
estimate  of  the  standard  deviation,  the  expected  half-width  is 
2.562//?  =  1.145  .  Finally,  assuming  the  spectral  method  with 

K  =  25  and  d  =  2  produces  an  unbiased  estimate  of  the  point 

2  1/2 

estimates'  standard  deviation,  (p(0;x  ) /mf  (x  )  '  =  1//?  . 

P  P 

Then  the  confidence  intervals  expected  half-width  is 
2.214//?  =  .99,  since  the  effective  number  of  degrees  of  freedom 
for  the  variance  estimate  was  shown  in  Table  1  of  Heidelberger 
and  Welch  (1981)  to  be  7. 

Thus,  given  the  assumptions,  the  spectral  confidence 
intervals  (with  d  =  2)  will  on  the  average  be  narrower  than  both 
the  agq  and  ngq  intervals  and  the  agq  intervals  will  be  slightly 
narrower  than  the  ngq  intervals.  If  a  cubic  polynomial  is  used 
in  the  spectral  method  then  the  spectral  confidence  intervals  can 
be  shown  to  be  somewhat  wider  than  the  agq  and  ngq  intervals. 

These  relationships  are  empirically  verified  in  Section  4. 


*■* •  .  - '*  • 


26 


4.  EMPIRICAL  STUDIES 


All  the  quantile  estimation  methods  given  in  the  previous 
section  are  asymptotically  valid  under  the  ^-mixing  assumption. 
The  crux  of  the  matter  though  is  whether  they  work  with  sample 
sizes  required  for  the  usual  precisions  required  in  systems  simu¬ 
lations,  say  10%  to  5%  as  measured  by  ratio  of  the  confidence 
interval  half-width  to  the  point  estimate.  In  this  section  we 
use  empirical  sampling  (simulation)  to  study  the  bias  and  stand¬ 
ard.  deviations  of  the  three  quantile  point  estimators  and  the 
confidence  interval  coverages  and  widths  of  the  three  quantile 
interval  estimators  described  in  Section  3.  The  tests  were  con¬ 
ducted  on  stationary  sequences  of  correlated,  exponentially  dis¬ 
tributed  random  variables  (NEAR (1)  and  GNEAR(l)  processes;  see 
Lawrance  and  Lewis  (1981a)  and  (1981b))  and  on  waiting  time  se¬ 
quences  in  heavily  congested  single  server  queues.  For  each 
process  the  0.50,  0.90,  0.99  and  0.999  quantiles  were  estimated. 

4.1.  Processes  Simulated 

The  NEAR ( 1 )  process  (Xn)  with  parameters  a  and  8 
is  defined  as  follows.  Let  {En>  be  a  sequence  of  i.i.d.  ex¬ 
ponentials  with  mean  1  and  let  {lnl  and  {Knl  be  mutua^y 
independent  sequences  of  random  variables  for  which 
P{Kn  =  1}  =  1  -  P(Kn  =  (l-a)6l  =  6,  where  6  =  ( 1— B ) / (1- (1-a) 8) , 
and  P{ln=8}=l-P{In=0}=a.  Set  X0  =  EQ  ;  then 
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is  a  Markovian  sequence  of  exponentially  distributed  random 
variables  with  mean  one  and  correlation  structure 

(4.2)  pk  =  corr(Xn,Xn+k)  =  (a6)k  =  pk  ,  k  >  0,  n  >  0 

The  GNEAR(l)  process  yields  a  process  with  exponential 
marginals  and  alternating  positive  and  negative  correlations. 
This  process  is  defined  by  XQ  =  Eg  and 

(4 . 3)  X  =  K  E  +  I  X'  ,  n  >  0  , 

n  n  n  n  n~l 

where  =  -  £n(l  -  exp(-  xn-1) )  •  For  the  GNEAR(l)  process 
the  lag  one  correlation  is 

(4.4)  =  (a6)r  , 

2 

where  r  =  1  -  IT  /6  =  -  0.6449  is  the  maximum  negative  cor¬ 
relation  attainable  in  a  bivariate  exponential  distribution. 
Efficient  algorithms  for  generating  the  NEAR ( 1 )  and  GNEAR(l) 
are  given  in  Lawrance  and  Lewis  (1981b). 

The  second  class  of  processes  we  considered  were 
(stationary)  waiting  time  sequences  in  M/G/l  queues,  where 
the  service  time  distribution  has  a  hyperexponential  distribu¬ 
tion.  Specifically,  let  Wn  denote  the  waiting  time  of  the 
nth  customer  and  let  {An  ,  n  ^  1}  and  {Sn  ,  n  0}  be  the 
i.i.d.  sequences  of  interarrival  times  and  service  times  re¬ 
spectively.  We  assume  that 
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1 


(4.5) 


P{ A  <  x>  =  1  -  exp (-  Xx)  , 
n  — 


x  >  0  ,  X  >  l 


(4.6)  p{Sn  i  x}  =  ~  exP(“  M^x)  }  +  n 2 { 1  -  exp(-  p2x)  )  , 


where  p^  and  are  greater  than  or  equal  to  zero,  and 

II i  +  n2  =  1  .  Let  p  =  (H1/p1  +  n2/p2)  1  and  let  p  =  X/p  <  1 
be  the  traffic  intensity.  Let  Wq  have  the  stationary  waiting 
time  distribution  (see  Kleinrock  (1975) ,  p,  205;  this  distribu¬ 
tion  is  a  probabilistic  mixture  of  an  atom  at  0  and  two  ex¬ 
ponentials)  ;  then  the  waiting  time  sequence  {Wn)  is  defined 
by 


(4.7) 


K„+l  '  <Wn  +  Sn  '  Vl> 


n  >  0 


where  x  =  max(x,0)  . 

We  considered  five  processes  in  testing  the  estimators. 
They  are: 


(i) 

NEAR ( 1 )  : 

a  = 

3 

=  0.95; 

(ii) 

NEAR  ( 1)  : 

a  = 

3 

=  0.995; 

( iii ) 

GNEAR(l)  : 

a  = 

3 

=  0.995; 

(iv) 

M/M/1 

X  = 

9, 

p^  =  10,  nx  =  1,  p  =  0.90; 

(v) 

M/G/l  : 

X  = 

9, 

p^  =  2 ,  p 2  =  18,  =  0.10 

(the  squared  coefficient  of  variation  of  the  service 
time  distribution  is  4.556). 
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The  correlation  structures  of  these  processes  range  from  moderate 
positive  correlation  (NEAR(l) ,  a  =  3  =  0.95,  Pk  ~  .9k)  to  extreme 
positive  correlation  (NEAR(l),  a  =  6  =  0.995,p  k  ~  .99k  ;  M/M/1 
and  M/G/l)  to  extreme  negative  correlation  (GNEAR(l) ) .  These 
processes  were  all  simulated  using  the  LLRANDOM  II  random  number 
generating  package  (see  Lewis  and  Uribe  (1981)). 


4.2.  Simulation  Results 

Tables  2-9  report  the  results  of  extensive  simulation 
studies  of  the  quantile  estimation  procedures  detailed  in  earlier 
Sections.  The  notation  used  in  the  tables  is  defined  as  follows. 
For  each  process,  quantile  level  p  ,  value  v  for  the  max- 
transform,  and  run  length  N  ,  R  i.i.d.  replications  were  per¬ 
formed  (R  *  200  for  all  runs,  except  runs  in  which  p  =  0.999, 
and  runs  in  which  p  =  .99  with  large  values  of  m) .  However,  the 
same  random  number  seeds  were  used  for  all  processes,  p,  v  and 
N  so  that  different  entries  in  the  tables  are  not  independent. 

A 


Let  Yq(r),  ngq(xp 

,  r)  and 

agq(xp  , 

r)  denote  the 

realiza- 

tions  on  the 

rth 

/\ 

replication  of  y^ 

,  ngq(xp) 

and 

agg(xp) 

respectively. 

Let 

yq  »  ngq 

and  agq 

denote 

the  averages 

over  the  R 

replications  of 

yq(r)  , 

ngq  (xp,r) 

and 

agq(xp,r) 

respectively. 

For 

example 

(4.8) 

R 

yq  =  I 

H  r=l 

yq(r)/R 

• 

Let  sd(yg)  , 

sd(ngq)  and  sd(agq)  denote  the 

sample 

i  standard 

deviations  of 

,  ngq  and 

agq  respectively 

.  For 

example 
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(4.9)  sd2(y  )  =  [  (y  (r)  _  y_)2/R(R-l)  . 

q  r=l  q  q 

These  empirical  sampling,  point  estimates  and  their  re-  | 

spective  estimated  standard  deviations  may  be  used  to  study  the  !  » 

bias  in  the  quantile  estimates  derived  in  this  paper.  For 

example,  from  Table  2,  an  approximate  90%  confidence  interval  j 

for  E(y  )  in  the  NEAR ( 1)  process  with  a  =  6  =  0.95,  p  =  0.99,  j 

q  i 

N  =  69000  and  v  =  69  is  y  +  1.645  sd(y  )  =  4.599  +  1.645  x  .009  .  ' 

q  —  q  — 

The  true  value  x  =  4.605  lies  within  this  90%  confidence  ! 

interval,  which  is  (4.584,  4.614). 

(i)  Table  2;  point  estimates  of  quantiles 

Table  2  illustrates  the  effect  of  the  maximum  transfor-  j 

mation  on  y  and  its  standard  deviation,  sd(y  )  .  For  each 

q  q  j 

process,  run  length  N  and  p  =  0.90  and  0.99,  y  and  sd(y  )  ' 

q  q  i 

were  computed  using  v  =  1  (no  max  transform)  and  v  chosen 

so  that  pv  ~  0.50.  The  run  lengths  N  were  chosen  so  that 

m  =  N/V  =  1000  for  the  NEAR ( 1 )  and  GNEAR(l)  processes  and  i 

m  =  N/v  =  2000  for  the  M/M/1  and  M/G/l  queues.  These  yield 

conservative  values  of  m.  The  only  case  in  which  the  max- 

transform  introduces  apparent  bias  is  the  M/M/1  queue  with 

p  =  0.90,  v  =  7,  N  =  14000;  in  this  case  y  differs  from  x 

P  P 

by  2.45  times  the  estimated  standard  deviation.  Since  this  is 
the  maximum  deviation  of  20  experiments,  it  is  probably  not  -j 

significant.  •  < 

The  column  labeled  "Actual  sd  inflation"  is  (with  p 
and  N  fixed)  the  ratio  sd(y  |pv  ~  0.50)/sd(y  |v=l)  .  This 

q  q 

ratio  measures  the  increase  in  standard  deviation  in  dependent 
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sequences  due  to  the  max- transform.  For  example,  for  the  NEAR ( 1) 
with  a  =  0  =  0.995  ,  p  =  0.99  and  N  =  69000  this  ratio  is 
.0090/. 0087  =  1.03  . 

The  column  labeled  "Theoretical  sd  inflation  for  i.i.d. 
samples"  is  the  inflation  in  standard  deviation  due  to  the  max- 
transform  for  a  sequence  if  i.i.d.  random  variables;  see 
equation  (2.8).  Recall  that  the  max-transform  changes  the  cor¬ 
relation  structure  of  the  process  so  that  its  effect  on  the 
standard  deviation  in  dependent  sequences  may  differ  substantially 
from  its  effect  in  independent  sequences.  In  most  cases  the 
actual  inflation  in  variance  is  very  slight  and  in  several  cases 
a  small  variance  reduction  is  achieved.  Only  for  the  M/M/1 
queue  is  the  inflation  greater  than  what  it  would  be  for  an  i.i.d. 
sequence,  and  even  in  this  case  it  is,  considering  the  storage 
and  computational  savings  of  the  max-transform, an  acceptable 
1.35  . 

From  Table  2  we  conclude  that,  for  sensibly  large  values 
of  m  ,  the  max-transform  introduces  very  little,  if  any,  bias 
and  that  the  inflation  in  variance  is  modest. 

(ii)  Tables  3-8;  confidence  interval  estimates 

Tables  3-8,  in  addition  to  reporting  the  estimated  means 
and  standard  deviations  of  the  quantile  point  estimates, 
compile  the  results  of  confidence  interval  coverage  studies  for 
the  three  quantile  confidence  interval  procedures.  For  each 
process,  p  ,  v  ,  N  and  replication  number,  confidence  inter¬ 
vals  for  Xp  were  generated  using  the  spectral  method,  nested 
group  quantiles,  and  averaged  group  quantiles  as  described  in 
Section  3. 
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For  these  studies,  G  ,  the  number  of  groups  was  set  at 
5  .  For  this  value  of  G  ,  under  the  assumptions  of  Section  3 
the  ngq  confidence  interval  has  nominal  coverage  0.9375,  i.e. 

A  A 

P^yq(l)  -  Xp  <  yq(5)^  =  0.9375.  For  comparison  purposes  we  also 
formed  93.75%  confidence  intervals  using  the  spectral  and  agq 
confidence  intervals.  Let  CIr(y  )  ,  cir (ngq)  and  CIr(agq) 
denote  the  confidence  interval  on  the  rth  replication  using 
the  spectral,  ngq  and  agq  methods  respectively  and  let 

/A 

|CIr(y^) |  ,  Ici^fngq) |  and  Ici^fagq) |  denote  the  widths  of 

these  intervals.  For  each  confidence  interval  procedure 
Tables  3-8  report  the  estimated  average  confidence  interval 
relative  half-width  (labeled  hw/x^)  .  For  example,  for  the 
spectral  confidence  intervals 

R 

(4.10)  hw/x  =  (1/R)  l  | Cl  (y  ) | /2x  . 

P  i=l  x  4  P 

These  tables  also  report  the  fraction  of  these  confidence 
intervals  which  actually  contain  xp  .  This  fraction  is  called 
an  (estimated)  93.75%  coverage  and  it  should  be  close  to 
0.9375  if  valid  confidence  intervals  are  being  formed.  A  one 
sided  confidence  interval  for  a  coverage,  based  on  the  normal 
approximation  to  the  binomial  distribution,  can  be  used  to  test 
the  hypothesis  that  the  actual  coverage  is  less  than  0.9375  . 

For  R  =  200  (Tables  3, 4, 5, 6  and  part  of  8)  estimated  coverages 
less  than  or  equal  to  0.91  are  significantly  low  at  the  0.90 
level  while  for  R  =  100  (Table  7  and  part  of  8)  estimated 
coverages  less  than  or  equal  to  0.89  are  significantly  low  at 
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the  0.90  level.  However,  in  practice,  the  importance  of  cover¬ 
age  in  judging  the  quality  of  a  confidence  interval  procedure 
depends  very  much  on  the  accuracy  (as  measured  by  relative  half¬ 
width)  of  the  confidence  intervals.  Thus  low  coverage  is  not  a 
serious  drawback  if  the  procedure  generates  confidence  intervals 
which  are  very  wide  (very  inaccurate)  and  therefore  provide  very 

limited  information  about  the  quantity  being  estimated.  In  fact, 

• 

the  most  useful  information  such  a  confidence  interval  often 
provides  is  that  the  run  length  is  too  short.  As  the  relative 
half-width  decreases,  the  importance  of  coverage  for  the  va¬ 
lidity  of  a  procedure  increases.  We  note  that  it  would  also  be 
possible  to  plot  the  coverage  functions  for  the  spectral  and 
agq  confidence  intervals  (Schruben  (1980));  however,  reporting 
the  individual  93.75%  coverages  is  representative  of  the  cover¬ 
age  function  and  more  compact. 

For  the  spectral  method  and  the  M/M/1  and  M/G/l  queues 
the  coverages  using  both  degrees  d  =  2  and  d  =  3  are  given 
whereas  for  the  NEAR ( 1 )  and  GNEAR(l)  processes  only  d  =  2  is 
given.  For  example,  from  Table  4,  the  coverages  for  xp  in  the 
M/M/1  queue  with  p  =  0.90,  v  =  1  and  N  =  14000  for  the 
spectral  method  using  degrees  d  =  2  and  3  ,  nested  group 
quantiles  and  averaged  group  quantiles  are  .905,  .910,  .880 
and  .870  respectively.  The  average  relative  half-widths  are 
.549,  .570,  .429  and  .403  respectively.  Thus  to  increase  the 
precision  from  its  current  value  of  about  50%  to  a  desired  value 
of  10%,  we  estimate  that  the  sample  size  would  have  to  be  in¬ 
creased  by  a  factor  of  25  to  a  total  of  350,000  observations. 


34 


It  would  clearly  be  impractical  to  store  and  sort  all  of  these 
observations.  Thus  to  increase  the  precision  to  an  acceptable 
level  requires  the  max- transform  to  reduce  the  sample  to  a 
manageable  size. 

Notice  that  when  the  coverages  of  all  methods  are  valid 
(for  example  in  Table  5  for  all  run  lengths  of  the  NEAR ( 1)  and 
GNEAR(l)  processes  and  large  run  lengths  of  the  queueing  processes) 
that  the  relative  widths  of  the  confidence  intervals  are  generally 
as  predicted  in  Section  3,  namely  that 

hw(spectral,  d=2)  £  hw(agq)  £  hw(ngq)  <  hw(spectral,  d=3)  . 

Tables  3  and  4  list  results  of  experiments  for  estimat¬ 
ing  the  median  and  0.90  quantile  respectively  without  using  the 
max-transform  (v  =  1) .  The  run  lengths  used  were  small  enough 
to  accommodate  storing  and  sorting  the  full  sequences. 

In  Table  3  all  coverages  are  either  acceptable  (i.e.f 
not  significantly  less  than  0.9375)  or  nearly  acceptable  (i.e., 
wide  confidence  intervals  with  coverage  close  to,  but  signifi¬ 
cantly  less  than  0.9375)  with  the  exception  the  GNEAR(l)  process 
using  the  spectral  method.  The  reason  for  the  low  coverage  in 
this  case  is  that  the  sequence  of  indicator  functions  {I^fx  5 ) } 
is  very  nearly  a  deterministic  sequence  of  alternating  zeros 
and  ones  due  to  the  strong  alternating  negative  and  positive 
correlation  structure  of  the  GNEAR(l)  process.  This  phenomenon 
appears  to  be  peculiar  to  the  median  (the  coverage  for  the  0.90 
quantile  using  the  spectral  method  in  Table  4  is  acceptable)  but 
should  be  kept  in  mind  when  dealing  with  processes  having  this 
type  of  correlation  structure. 
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The  coverages  for  the  0.90  quantiles  in  Table  4  are  all 


acceptable  or  nearly  acceptable.  However,  agq  and  ngq  exhibit 
substantial  bias  in  the  M/G/l  queue.  In  both  Tables  3  and  4 
notice  the  very  large  mean  relative  half-widths  in  the  M/M/1 
and  M/G/l  queues,  again  implying  the  need  for  the  max-transform 
to  deal  with  the  much  larger  sample  sizes  required  for  estimates 
of  acceptable  precision. 

Tables  5,  6  and  7  compile  results  of  experiments  for 
estimating  the  0.90,  0.99  and  0.999  quantiles,  respectively, 
using  the  max-transform  with  v  chosen  so  that  pv  ~  0.50  . 

The  ngq  and  agq  confidence  interval  estimate  coverages 
are  generally  acceptable  throughout  these  tables.  However,  as 
in  the  case  of  v  =  1  ,  ngq  and  agq  show  that  there  is  sub¬ 
stantial  bias  in  the  ngq  and  agq  point  estimates  in  the 
M/G/l  queue  for  small  N.  For  N  =  224,000  this  bias  has  dis¬ 
appeared  even  though  the  precision  (about  20%)  is  still  relatively 
high.  The  estimates  x^  show  little,  if  any,  bias.  The  cover¬ 
ages  of  the  confidence  interval  estimates  in  the  M/G/l  queue 
are  high  despite  the  bias  due  to  the  extreme  width  of  the  confi¬ 
dence  intervals  at  the  small  run  lengths. 

For  the  spectral  method  the  coverages  for  the  NEAR(l) 
and  GNEAR(l)  processes  are  acceptable.  However,  the  small  sample 
coverages  for  the  queues,  particularly  the  M/G/l  queue,  are 
disappointing.  Even  using  a  cubic  polynomial  to  fit  the  log  of 
the  periodogram  does  not  yield  acceptable  coverages  for  the 
M/G/l  queue  when  N  is  small.  The  low  coverage  is  explained 
by  the  fact  that  the  spectrum  of  the  max-transformed  variables 
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{ Yk }  is  still  very  steep  near  f  =  0  (i.e.,  the  Y^s  are 

highly  correlated)  even  for  large  values  of  v  .  Since  N  =  mv  , 
if  v  is  large  then  m  is  small  and  the  least  squares  fits  to 
the  log  of  the  periodogram  is  done  over  a  relatively  large  fre¬ 
quency  range.  If  log  (p(f  ;  x  )  )  is  very  steep  then  it  may  not 
be  well  approximated  by  a  low  order  polynomial  over  a  wide  fre¬ 
quency  range. 

Figure  1  illustrates  this  problem  in  the  case  of  the 
NEAR(l)  process  with  a  =  £  =  0.995  .  The  log  of  the  spectrum 
of  the  max-transformed  variables  (v  =  693)  is  seen  to  be  still 
quite  steep  near  f  =  0  .  Compare  this  to  the  almost  flat  log 
of  pn(f)  ,  the  spectrum  of  the  process  of  batch  means  X_.(k)  , 

where 


(4.11) 


Vk)  =  i 


kB 

l 

j  = (k-1) B+l 


These  batch  means  would  be  used  to  place  a  confidence  interval 
on  E  (X j )  (see  Heidelberger  and  Welch  (1981)). 

The  problem  can  be  avoided  by  making  v  smaller,  and 
thus  m  larger.  The  least  squares  fit  is  then  done  over  a 
narrower  frequency  range  and  the  coverage  improves.  This  effect 
is  seen  by  comparing  the  coverages  for  x  in  the  M/G/l 

queue  with  N  =  14000  and  v  =  1  (Table  4)  to  those  using  v  =  7 
and  N  =  14000  (Table  5) .  (Again  note  that  the  precision  is 
unacceptably  high).  It  is  also  seen  in  Table  8  which,  for  the 
M/M/1  and  M/G/l  queues,  compares  the  coverages  for  x  „„ 
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log  spectrum 


V. 


J 

f. 


Frequency,  f 

Figure  1.  Effect  of  the  maximum  transform  and  batching 
on  the  log-spectrum.  Figure  la  gives  the  log 
of  the  spectrum  of  a  NEAR(l)  process  with 
a  =  6  =  0.995.  Figure  lb  shows  the  log  of 
the  spectrum  after  the  maximum  transformation 
with  v  =  693.  A  small  but  significant  effect 
can  be  seen  at  f  =  0.  Figure  lc  shows  the 
log  of  the  spectrum  after  batching,  with  batch 
sizes  of  v  =  693.  Batching  reduces  the  initial 
point  of  the  spectrum  much  more  than  the  maximum 
transformation. 
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log  spectrum  of  max  transformed  variables  (v=693) 


Figure  1.  Effect  of  the  maximum  transform  and  batching 
on  the  log-spectrum.  Figure  la  gives  the  log 
of  the  spectrum  of  a  NEAR(l)  process  with 
a  =  3  =  0.995.  Figure  lb  shows  the  log  of 
the  spectrum  after  the  maximum  transformation 
with  v  =  693.  A  small  but  significant  effect 
can  be  seen  at  f  =  0.  Figure  lc  shows  the 
log  of  the  spectrum  after  batching,  with  batch 
sizes  of  v  =  633.  Batching  reduces  the  initial 
point  of  the  spectrum  much  more  than  the  maximum 
transformation. 
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log  spectrum  of  batch  means.  Batch  size  =  693 


Figure  1.  Effect  of  the  maximum  transform  and  batching 
on  the  log-spectrum.  Figure  la  gives  the  log 
of  the  spectrum  of  a  NEAR(l)  process  with 
a  =  6  =  0.995.  Figure  lb  shows  the  log  of 
the  spectrum  after  the  maximum  transformation 
with  v  =  693.  A  small  but  significant  effect 
can  be  seen  at  f  =  0.  Figure  1c  shows  the 
log  of  the  spectrum  after  batching,  with  batch 
sizes  of  v  =  693.  Batching  reduces  the  initial 
point  of  the  spectrum  much  more  than  the  maximum 
transformation. 


40 


I 


I 

! 


i 


using  v  =  69  to  those  using  v  =  17  and  v  =  8  ,  and  the 
coverages  for  x  ggg  using  v  =  693  to  those  using  v  =  173 
and  v  =  86  (the  slight  differences  in  N  are  to  make  m 
highly  composite  which  increases  the  computational  efficiency  of 
the  Discrete  Fourier  Transform  used  to  calculate  the  periodogram) . 
Thus  the  full  storage  and  computational  savings  of  the  max- 
transform  (i.e.  moving  the  quantile  back  to  the  median)  cannot 
be  realized  in  highly  correlated  sequences  without  sacrificing 
confidence  interval  coverage.  Significant  savings  can  be  realized 
however.  For  example  the  max-transform  with  v  =  173  reduces 
the  1,384,000  observations  to  a  sample  size  of  just  8000  for 
estimating  x  ggg  in  the  M/G/l  queue. 

In  Table  8  two  simulations  are  also  shown  with 
N  =  5,536,000  and  v  =  173  for  both  the  M/M/1  and  M/G/l 
queues  for  the  extreme  case  of  the  .999  quantiles.  This  sample 
size  is  large  enough  to  obtain  almost  10%  precision.  To  within 
the  precision  of  the  empirical  sampling  experiment  (R  =  100) ,  it 
should  be  noted  that  all  point  and  confidence  interval  estimates 
are  valid.  Note  that  the  maximum  transformation  reduces  the 
sample  size  from  5,536,000  to  32,00011 
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5.  SUMMARY  AND  CONCLUSIONS 

In  this  paper  we  have  investigated  point  and  interval 
estimates  for  a  single  quantile  in  a  fixed  length  sequence  of 
dependent  observations.  For  large  sample  sizes  it  becomes  im¬ 
practical  to  store  and  sort  the  entire  sequence.  For  extreme 
quantiles,  these  limitations  can  be  overcome  by  using  the  maxi¬ 
mum  transformation  which  requires  storing  only  a  sequence  of 
maxima.  This  sequence  is  defined  by  laying  the  data  out  into  a 
v  by  m  array  and  storing  only  the  maximum  element  in  each 
column.  Storage  requirements  are  thus  reduced  by  a  factor  of  v  . 
Observations  at  lag  m  are  assumed  to  be  independent  and  the 
max-transform  changes  the  problem  of  estimating  the  p  quantile 
of  the  original  sequence  into  one  of  estimating  the  pV  quantile 
of  the  max- transformed  sequence. 

Three  confidence  interval  methods  which  can  exploit  the 
sample  size  reduction  produced  by  the  max-transformation  were 
described  and  tested;  the  spectral  method  and  two  extensions  to 
the  method  of  batch  means  called  nested  group  quantiles  (ngq) 
and  averaged  group  quantiles  (agq) .  The  three  methods  produced 
valid  confidence  intervals  if  the  sample  sizes  were  large  enough 
to  produce  10%  precision  in  the  estimates. 

Problems  which  were  not  addressed  in  this  paper  and  are 
the  subject  of  ongoing  research  involve  questions  of  multiple 
quantile  estimation  and  the  incorporation  of  these  fixed  run 
length  procedures  into  sequential  procedures  for,  say,  simulation 
run  length  control.  Among  the  issues  involved  here  are: 
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(i)  Testing  the  adequacy  of  the  spacing  parameter  m 
so  that  and  xk+m  are  approximately  independent. 

(ii)  Organizing  the  data  so* that  m  can  be  increased 
as  the  sample  size  increases. 

(iii)  Organizing  the  data  so  that  multiple  quantiles  can 
be  efficiently  estimated,  for  example,  estimating  a  0.50  and 
0.999  quantile  for  the  same  set  of  observations  (this  probably 
involves  the  use  of  a  max-min  transformation) . 

(iv)  Determining  the  sensitivity  of  these  schemes  to 
an  initial  transient  which  is  often  encountered  in  simulation. 
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Table  2.  Comparison  of  point  estimates  and  standard  deviations  of  point  estimates 

using  full  sample  order  statistic  and  max-transformed  sample  order  statistic 
(R  -  200  replications).  Here  q  -  pv  .  The  fact  that  the  actual  sd  inflation 
is  generally  less  than  the  theoretical  sd  inflation  for  i.i.d.  sequences  is 
because  of  the  decrease  in  correlation  in  the  series  induced  by  the  maximum 
transformation. 


— 

Process 

P 

X 

p 

N 

1 

sd(y^) 

Actual 

sd 

inflation 

Theoretical  sd 
inflation  for 
i.l.d.  sequences 

NEAR(l) 

2.303 

7p00 

i 

a- 6-0.  95 

2.303 

7P00 

7 

2.320 

B 

1.02 

1.18 

4.605 

69P00 

1 

4.610 

.0087 

4.605 

69 

4.599 

.0090 

1.03 

1.20 

— 

NEAR(l) 

0.90 

2.303 

7,000 

1 

2.344 

19 

a- 6*0.995 

0.90 

2.303 

7p00 

7 

2.340 

EH 

1.10 

1.18 

0.99 

4.605 

69P00 

1 

4.540 

.0275 

0.99 

4.605 

69 

4.607 

.0299 

1.09 

1.20 

GNEAR(l) 

0.90 

2.303 

7,000 

■ 

2.282 

a-6*0. 995 

0.90 

2.303 

7P00 

B 

2.291 

0.92 

1.18 

negative 

M 

correlation 

0.99 

4.605 

69P00 

■9 

4.605 

.0046 

0.99 

4.605 

69000 

69 

4.600 

.0050 

1.09 

1.20 

M/M/1 

0.90 

2.197 

14P00 

M 

2.222 

rib 

P-0.90 

0.90 

2.197 

2.289 

mm 

1.35 

1.18 

0.99 

4.500 

138P00 

■ 

4.514 

.0368 

0.99 

4.500 

139000 

69 

4.558 

.0455 

1.24 

1.20 

M/G/l 

0.90 

6.311 

14000 

1 

.1682 

P-0.90 

0.90 

6.311 

14P00 

7 

.1759 

1.05 

1.18 

0.99 

13.130 

B 

1 

13.182 

0.99 

13.130 

B 

69 

13.398 

0.95 

1.20 

V, 


Table  3.  Point  estimates,  standard  deviations  of  point  estimates,  93.75% 
coverages,  and  average  confidence  interval  half-width/x  for  the 

p 

three  methods.  Here  p  =  0.50,  v  -  1,  G  =  5,  and  R,  the  number 
of  replications,  is  200. 


Process 

X 

p 

N 

Spectral  Method 

Nested  Group 
Quantiles 

Averaged  Group 
Quantiles 

X 

P 

sd(x  ) 

P 

degree  of  polynomial 

3  2 

coverage 

hw/x 

P 

coverage 

hw/x 

P 

coverage 

hw/x 

P 

coverage 

hw/x 

P 

ngq 

sd(ngq) 

ngq 

sd(agq) 

NEAR(l) 
a- 6-0. 95 

0.693 

7p00 

. 

.696 

(.003) 

.935 

(.149) 

.700 

(.004) 

.935 

(.180) 

mm 

.945 

(.167) 

I 

NEAR(l) 
a- 6*0. 995 

0.693 

7,000 

.865 

(.517) 

.744 

(.015) 

.905 

(.598) 

.809 

(.013) 

.910 

(.565) 

GNEAR(l) 
ct- 6-0.  995 
(negative 
correlation] 

0.693 

7,000 

.692 

(.001) 

_  ... 

.795 

(.006) 

.693 

(.002) 

.945 

(.112) 

.985 

(.102) 

M/M/1 

P-0.90 

0.588 

14p00 

.601 

(.006) 

.950 

(.452) 

.950 

(.352) 

.598 

(.007) 

.945 

(.443) 

.640 

(.008) 

.940 

(.413) 

M/G/l 
p-0. 90 

1.5A5 

14p00 

1.607 

(.033) 

.950 

(.878) 

.940 

(.647) 

1.562 

(.036) 

.930 

(1.013) 

1.943 

(.055) 

.940 

(.961) 

Table  4.  Point  estimates,  standard  deviations  of  point  estimates,  93.75% 
coverages,  and  average  confidence  interval  half-width/x  for  the 
three  methods,  no  maximum  transform.  Here  p  -  0.90,  v  »  1, 

G  -  5  and  R,  the  number  of  replications,  is  200. 


Process 

X 

p 

N 

NEAR(l) 
a- 6-0. 95 

2.303 

7,000 

NEAR(l) 
a- 6-0. 995 

2.303 

7P00 

GNEAR(l) 
a- 6-0. 995 
(negative 
correlation) 

2.303 

7,000 

M/M/1 

p-0.90 

2.197 

upoo 

M/G/l 

p-0.90 

6.311 

14p00 

Spectral  Method 

degree  of  polynomial 

3  2 

X 

p 

coverage 

coverage 

sd(x  ) 

hw/x 

hw/x 

P 

P 

P 

(.010) 


.910 

(.570) 


.890 

(1.069) 


coverage 


coverage 


.930 

(.129) 


.935 

(.468) 


.925 

(.141) 


.905 

(.549) 


8d(ngq)  hw/Xjj  sd(agq)  hw/x^ 


2.296 

(.010) 


2.282 

(.013) 


.955 

(.181) 


.880 

(.429) 


5.266 

(.090) 


.940 

(.435) 


2.229 

(.027) 

.870 

(.403) 

5.843 

(.099) 

.815 

(.466) 

V 
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Table  5.  Point  estimates,  standard  deviations  of  point  estimates,  93.75% 
coverages,  and  average  confidence  interval  half-width/x  for  the 
three  methods.  As  in  Table  4,  p  =  .90,  but  here  v  =  7;P  thus  pv 
Also  G  =  5  and  R,  the  number  of  replications,  is  200. 


0.50  . 


Process 

X 

p 

N 

Spectral  Method 

Nested  Group 
Quantiles 

Averaged  Group 
Quantiles 

sd(yq) 

degree  of 

3 

polynomial 

2 

coverage 

hw/x 

P 

1 

coverage 

hw/x 

P 

coverage 

hw/x 

P 

coverage 

hw/x 

P 

ngq 

sd(ngq) 

ngq 

sd(agq) 

NEAR(l) 

2.303 

7,000 

2.  320 

.945 

2.310 

.940 

2.328 

.940 

a* 0*0. 95 

(.011) 

(.166) 

(.025) 

(.179) 

(.022) 

(.164) 

28,000 

2.304 

.940 

2.303 

.920 

2.305 

.915 

(.006) 

(.076) 

(.007) 

(.089) 

(.006) 

(.083) 

112,000 

2.300 

.930 

2.296 

.920 

2.300 

.935 

(.003) 

(.037) 

(.003) 

(.044) 

(.000) 

(.040) 

NEAR(l) 

2.303 

7,000 

2.  340 

.800 

2.253 

.940 

2.414 

.945 

a=0=O. 995 

(.036) 

(.342) 

(.038) 

(.  629) 

(.033) 

(.488) 

28,000 

2.311 

.935 

2.297 

.960 

2.350 

.940 

(.017) 

(.246) 

(.020) 

(.279) 

(.017) 

(.  200)  > 

112,000 

2.289 

2.301 

.945 

2.309 

.945 

(.008) 

(.127) 

(.010) 

(.142) 

(.008) 

(.  232)  ; 

GNEAR(l) 

2.303 

7,000 

2.291 

.960 

2.316 

.935 

2.287 

.940 

a- 0-0. 995 

(.012) 

(.199) 

(.013) 

(.292) 

(.022) 

(.278)  ! 

28,000 

2.306 

.940 

2,312 

.940 

2.299 

.935  ! 

(negative 

(.006) 

(.078) 

(.007) 

(.097) 

(.006) 

(.090)  j 

correlation) 

112,000 

2.304 

.925 

2.306 

.925 

2.304 

.935 

_ _ 

(.003) 

(.040) 

(.00J) 

(.047) 

(.003) 

(.043) 

M/M/1 

2.197 

2.289 

.940 

.815 

.955 

2.331 

.930 

P-0.90 

(.037) 

(.300) 

1.306) 

(.033) 

(.461) 

(.002) 

(.422) 

2.231 

.975 

.950 

2.220 

.945 

2.316 

.935 

(.016) 

(.326) 

(.240) 

( . 020) 

(.290) 

(.020) 

(.274) 

224,000 

2.206 

.965 

.950 

2.204 

.955 

2.227 

.945 

(.008) 

(.142) 

(.116) 

(.00$) 

.  . 

(. 203) 

(.008) 

(.227) 

M/G/l 

6.311 

14,000 

6.572 

.695 

.525 

5.657 

.865 

6.112 

.860 

c-0.90 

(.176) 

(.399) 

(.219) 

(.123) 

(.500) 

(.220) 

(.472) 

56,000 

6.373 

.950 

.870 

6.184 

.940 

6.680 

.935 

(.079) 

(.471) 

(.300) 

(.075) 

(.444) 

(.086) 

(.417) 

224,000 

6.388 

.970 

.955 

6.367 

.965 

6.559 

.960 

(.042) 

(.260) 

(.  209) 

(.047) 

(.243) 

(.047) 

(.227) 

i 


i 


’3 


* 


4 


*4 

% 

$ 


i 


II 
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Table  6.  Point  estimates,  standard  deviations  of  point  estimates,  93.75% 
coverages,  and  average  confidence  interval  half-width/x  for  the 
three  methods.  Here  p  «  0.99,  v  -  69,  so  that  ^ 

9  “  Pv  £  0.50  .  Also  G  -  5  and  the  number  of  replications,  R, 
is  200.  However  R  -  100  if  the  run  is  marked  t  . 


Process 

1 

N 

Spectral  Method 

Nested  Group 
Quantiles 

I 

Averaged  Group 
Quantiles 

*q 

sd(yq) 

degree  of  polynomial 

3  2 

coverage 

hw/x 

P 

coverage 

hw/x 

P 

coverage 

hw/x 

P 

ngq 

sd(ngq) 

ngq 

sd(agq) 

coverage 

hw/x 

P 

NEAR(l) 

4.605 

69,000 

4.599 

.960 

4.606 

.920 

4.627 

.915 

a- 6-0. 95 

(.  009) 

(.  062) 

(.011) 

(.073) 

(.009) 

(.002) 

138,000 

4.608 

.945 

4.604 

.950 

4.608 

.950 

{.006) 

(.042) 

(.007) 

(.054) 

(.006) 

(.050) 

276,000 

4.604 

.930 

4.603 

.945 

4.608 

.950 

{.004) 

(.021) 

(.005) 

(.027) 

(.005) 

(.034) 

NEAR(l) 

4.605 

69,000 

4.607 

.890 

4.606 

.965 

4.692 

.965 

{.020) 

(.180) 

(.032) 

(.248) 

(.029) 

(.229) 

138,000 

4.653 

.930 

4.653 

.905 

4.685 

.925 

{.022) 

(.140) 

(.027) 

(.172) 

(.024) 

(.160) 

276,000 

4.620 

.940 

4.613 

.945 

4.629 

.935 

{■014) 

(.099) 

(.016) 

(.111) 

(.014) 

(.702) 

GNEAR(l) 

4.605 

69,000 

4.600 

-  — 

.955 

4.596 

.950 

4.594 

.940 

a- 6-0. 995 

{.005) 

(.025) 

(.006) 

(.040) 

(.005) 

(.032) 

138,000 

4.605 

.925 

4.601 

.910 

(negative 

{.004) 

(.023) 

(.005) 

(.020) 

(.004) 

(.020) 

correlation) 

276,000 

.950 

4.604 

.930 

4.605 

.935 

{.002) 

(.017) 

(.003) 

(.020) 

(.003) 

(.019) 

M/M/1 

4.500 

138,000 

4.558 

.865 

.730 

4.351 

.920 

4.538 

.905 

{.045) 

(.262) 

(.153) 

(.035) 

(.267) 

(.035) 

(.250) 

276,000 

4.513 

.955 

.895 

4.472 

.930 

4.598 

.930 

{.027) 

(.226) 

(.148) 

(.032) 

(.222) 

(.020) 

(.202) 

552,000 

4.525 

.955 

.945 

4.511 

.945 

4.608 

.935 

{.020) 

(.190) 

(.132) 

(.024) 

(.178) 

(.023) 

(.166) 

t 

2,208,000 

4.481 

.980 

.960 

4.495 

.960 

4.510 

.950 

{.012) 

(.084) 

(.067) 

(.015) 

(.076) 

(.0.Z3) 

(.020) 

M/G/l 

13.13C 

138,000 

13.398 

.675 

.490 

11.713 

mmfm 

12.174 

.815 

{.198) 

(.242) 

(.129) 

(.125) 

(.119) 

(.220) 

276,000 

13.359 

.845 

.635 

12.377 

12.934 

.900 

{.154) 

(.270) 

(.155) 

(.107) 

(.287) 

(.106) 

(.202) 

552,000 

13.271 

.935 

.815 

12.927 

.910 

13.480 

.915 

{.106) 

(.276) 

(.164) 

(.103) 

(.205) 

(.101) 

(.242) 

t 

2,208,000 

[13.102 

.960 

.930 

13.010 

.920 

13.241 

.900 

{.062) 

(.163) 

(.123) 

(.075) 

(.122) 

(.  000) 

(.125) 

.  •  .  -v  ,* 


Table  7.  Point  estimates,  standard  deviations  of  point  estimates,  93.75% 
coverages,  and  average  confidence  interval  half-width/x  for  the 
three  methods.  Here  p  =  0.999,  v  =  693,  so  that  ^ 

q  =  pv  ~  0.50  .  Also  G  -  5  and  the  number  of  replications,  R, 
is  100. ~ 


Spectral  Method 


NEAR(l)  6.908  693000 

a=6=0. 95 


NEAR(l)  6.908  693000 

a- 6=0. 995 


GNEAR(l)  6.908  693000  6.9 

a=6=0. 995  (.0 

(negative 
correlation) 


degree  of  polynomial 

3 _  2  _  _ 

coverage  coverage  ngq  coverage  ngq  coverage 

hw/Xp  hw/x^  sd(ngq)  hw/x^  sd(agq)  hw/x^ 


.940  6.908  .940  6.912  .940 

(.034)  (.012)  (.038)  (.010)  (.035) 


.930 

(.016) 
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Table  8.  M/M/1  queue.  Point  estimates,  standard  deviations  of  point 

estimates,  93.75%  coverages,  and  average  confidence  interval 
half-widths/x  comparing  q  =  pv  ~  0.50,  q  =  pv  ~  0.84  and 
9  “  pv  ~  0.92  for  different  sample  sizes.  Here  p  «  0.99  and 
p  •  0.999  and  G  =  5.  The  number  of  replications  is  R  =  200 
for  runs  marked  *  ,  and  R  =  100  for  runs  marked  t  . 


Spectral  Method 

Nested  Group 
Quantiles 

Averaged  Group 
Quantiles 

_ 

degree  of  polynomial 

3  2 

P 

y9 

coverage 

coverage 

ngq 

coverage 

ngq 

coverage 

Process 

X 

p 

N 

sd(y^) 

hw/ x 

P 

hw/x 

P 

sd(ngq) 

hw/x 

P 

sd(agq) 

hw/x 

P 

M/M/1 

p-0.90 

* 

0.99 

4.500 

138,000 

69 

4.558 

(.045) 

.865 

(.262) 

.730 

(.153) 

4.351 

(.035) 

.920 

(.267) 

4.538 

(.035) 

.905 

(.250) 

* 

0.99 

4.500 

136,000 

17 

4.546 
(.  040) 

.925 

(.375) 

.880 

(.246) 

4.258 

(.029) 

.895 

(.240) 

4.409 

(.029) 

.890 

(.222) 

* 

128,000 

8 

4.596 

(.042) 

.900 

(.384) 

.870 

(.306) 

4.317 

(.035) 

.895 

(.260) 

4.470 

(.037) 

.875 

(.243) 

t 

0.999 

6.800 

1,386,000 

693 

6.835 

(.058) 

CJE,  Bq 

1 

| 

.930 

(.172) 

6.790 

(.046) 

.940 

(.160) 

+ 

0.999 

6.800 

1,384,000 

173 

6.855 

(.059) 

.910 

(.284) 

.870 

(.178) 

6.578 

(.055) 

.930 

(.178) 

6.732 

(.050) 

.910 

(.165) 

t 

0.999 

6.800 

1,376,000 

86 

6.856 

(.056) 

.940 

(.232) 

.890 

(.181) 

6.657 

(.054) 

.900 

(.156) 

6.754 

(.049) 

.880 

(.144) 

t 

0.999 

6.800 

5,536,000 

173 

6.797 

(.024) 

.970 

(.112) 

.950 

(.088) 

6.743 

(.027) 

.930 

(.108) 

6.833 

(.028) 

.940 

(.100) 

M/G/l 
p=0. 90 

* 

0.99 

13.13C 

138,000 

69 

13.398 

(.198) 

.675 

(.243) 

.490 

(.129) 

11. 713 
(.125) 

.820 

(.288) 

12.174 

(.119) 

.815 

(.270) 

* 

0.99 

13.13C 

136,000 

17 

13.24c 

(.202) 

.850 

(.462) 

.725 

(.274) 

11.467 

(.125) 

.795 

(.294) 

12.011 

(.124) 

.755 
(.  277) 

* 

128,000 

8 

13.011 

(.200) 

.805 
(.  758) 

.730 

(.489) 

11.336 

(.126) 

.755 

(.286) 

11.762 

(.115) 

.720 

(.266) 

t 

0.999 

19.94 

1,386,000 

693 

20.693 

(.324) 

.  700 
(.161) 

18.589 

(.171) 

.890 

(.211) 

19.063 

(.157) 

.840 

(.198) 

t 

0.999 

19.94 

1,384,000 

173 

19.951 
(.  230) 

.890 

(.294) 

18.188 

(.159) 

.880 

(.198) 

18.638 

(.137) 

.880 

(.185) 

t 

0.999 

19.94 

1,376,000 

86 

19.844 

(.239) 

.820 

(.385) 

.780 

(.238) 

18.053 

(.176) 

.800 

(.198) 

18.621 

(.159) 

.770 

(.183) 

t 

0.999 

19.94 

5,536,000 

173 

20.152 

(.146) 

.950 

(.211) 

.920 

(.156) 

19.511 

(.138) 

.920 

(.151) 

19.979 

(.128) 

.950 

(.142) 
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